Italian preterms 16S and virome analysis

1 Sequencing depth

Both bacterial and viral samples have some very low-read samples - especially at 1 Days

2 Rarefaction curves

All samples have sufficient read depth, except for 37297 and 27281 - these sample should possibly be removed

## Warning: Could not find a column named OTU/ASV in otutable, using rownames as
## OTU ID's

## Warning: Could not find a column named OTU/ASV in taxonomy, using rownames as
## OTU ID's

## Warning: Could not find a column named OTU/ASV in otutable, using rownames as
## OTU ID's

## Warning: Could not find a column named OTU/ASV in taxonomy, using rownames as
## OTU ID's

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 5

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 57

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 8

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 2

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 9

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 2

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 29

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 11

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 9

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 9

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 4

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 2

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 9

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 99

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 9

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 2

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 2

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 5

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 3

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 2

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 9

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 7

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 5

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 3

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 13

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 2

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 18

## Warning in vegan::rarefy(abund[i, ], n): most observed count data have counts
## 1, but smallest count is 46

Do we have any positive or negative controls?

3 Viral taxonomy stats

total_reads <- sum(sample_sums(PSV))

#Realm-level
# no_realm <- subset_taxa(PSV,Realm=="Unclassified")
# 1-sum(sample_sums(no_realm))/total_reads

#Class
# no_class <- subset_taxa(PSV,Class=="Unclassified")
# 1-sum(sample_sums(no_class))/total_reads

#Order
no_order <- subset_taxa(PSV,Order=="Unclassified")
1-sum(sample_sums(no_order))/total_reads

## [1] 0.0689926

#Family
no_family <- subset_taxa(PSV,Family=="Unclassified")
1-sum(sample_sums(no_family))/total_reads

## [1] 0.005573007

#Genus
no_genus <- subset_taxa(PSV,Genus=="Unclassified")
1-sum(sample_sums(no_genus))/total_reads

## [1] 0

#Species
no_species <- subset_taxa(PSV,Species=="Unclassified")
1-sum(sample_sums(no_species))/total_reads

## [1] 0

#Caudoviricetes
caudovirales <- subset_taxa(PSV,Class=="Caudoviricetes")
sum(sample_sums(caudovirales))/total_reads

## [1] 0.994427

caudovirales_uncl <- subset_taxa(caudovirales,Genus=="Unclassified")
sum(sample_sums(caudovirales_uncl))/sum(sample_sums(caudovirales))

## [1] 1

sum(sample_sums(caudovirales_uncl))/total_reads

## [1] 0.994427

#Classified caudovirales
(sum(sample_sums(caudovirales))-sum(sample_sums(caudovirales_uncl)))/total_reads

## [1] 0

4 Barplots

4.1 Host predictions

5 Abundance heatmaps

5.1 Bacteria

5.2 Virome

5.3 Host virome

6 Upset plot

6.1 Bacteria

6.2 Virome

6.3 Virome host

7 Alpha diversity

7.1 Bacterial and viral shannon diversity

7.2 Bac vs viral alfa diversity

#Subset samples present in both
vir_ids <- PSV@sam_data$description %>% str_split_i("_",-1)
PSV@sam_data$id <- vir_ids
PSBR@sam_data$description <-PSBR@sam_data$description %>% str_split_i("_",-1)
PSVR@sam_data$description <-PSVR@sam_data$description %>% str_split_i("_",-1)

PSB.R.both <- subset_samples(PSBR,description %in% vir_ids)
#PSV.CSS.both <- subset_samples(PSV.CSS,vir_ids %in% PSB.CSS.both@sam_data$description)
PSV.R.both <- subset_samples(PSVR,vir_ids %in% PSB.R.both@sam_data$description)

shannon_bac = estimate_richness(PSB.R.both, measures = "Shannon")
shannon_bac <- cbind(shannon_bac, Days = sample_data(PSB.R.both)$Days) %>% dplyr::rename(Bac_Shannon = Shannon)

shannon_vir = estimate_richness(PSV.R.both, measures = "Shannon")
shannon_vir <- cbind(shannon_vir, Days = sample_data(PSV.R.both)$Days) %>% dplyr::rename(Vir_Shannon = Shannon)

shannon_bac_vir <- cbind(shannon_bac,Vir_Shannon=shannon_vir$Vir_Shannon)

ggplot(shannon_bac_vir,aes(x=Bac_Shannon,y=Vir_Shannon
                           ,color=Days
                           ))+
  geom_point() 

  # gpmisc::stat_poly_line() +
  # ggpmisc::stat_poly_eq(use_label(c("eq", "R2")))"

8 Beta diversity

8.1 Days

NB! This is including some bacterial day 0 samples >100 reads!

9 Within vs between twin pairs bray curtis dissimilarity

9.1 Individual time points

9.1.1 Day 1

9.1.2 Day 15

9.1.3 Day 30

9.1.4 Day 90

10 Effect size

10.1 Bacteria

## [1] 137

10.1.1 Db anova - constrained by time point

10.2 Virome

## [1] 576

10.2.1 Db anova - constrained for time point

10.3 Host virome

## [1] 305

10.3.1 Db anova - constrained for time point

10.4 Host virome - Genus

## [1] 305

##  [1] "Twin_pairs"                         "Days"                              
##  [3] "Gestational_age"                    "Birth_weight"                      
##  [5] "Small_for_gestational_age"          "Antenatal_steroids"                
##  [7] "Delivery_mode"                      "Intrapartum_antibiotic_prophylaxis"
##  [9] "Late_onset_sepsis"                  "Necrotizing_enterocolitis"         
## [11] "Patent_ductus_arteriosus"           "Bronchopulmonary_dysplasia"        
## [13] "Antibiotics_at_birth"               "Total_days_antibiotics_therapy"

10.4.1 Db anova - constrained for time point

10.5 Combined effect size plots

10.5.1 With viral hosts

10.5.2 Bacteria vs vira only

11 Rabuplot

11.1 Time

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

PSV.host@sam_data$Days <- factor(PSV.host@sam_data$Days, levels = c("90","60","30","15","1"))
PSB@sam_data$Days <- factor(PSB@sam_data$Days, levels = c("90","60","30","15","1"))

rabu.vir <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Days", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:4]) +
  ggtitle("Viral hosts")

rabu.bac <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Days", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:4]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria")

ggarrange(rabu.bac,rabu.vir, common.legend = TRUE, legend = "bottom")

11.2 Bronchopulmonary_dysplasia

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.bpd <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Bronchopulmonary_dysplasia", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:2]) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

rabu.bac.bpd <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Bronchopulmonary_dysplasia", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:2]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

ggarrange(rabu.bac.bpd,rabu.vir.bpd, common.legend = TRUE, legend = "bottom")

11.3 Intrapartum antibiotics

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.tdab <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Intrapartum_antibiotic_prophylaxis", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

rabu.bac.tdab <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Intrapartum_antibiotic_prophylaxis", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

ggarrange(rabu.bac.tdab,rabu.vir.tdab, common.legend = TRUE, legend = "bottom")

11.4 Total days antibiotics treatment

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.iap <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Total_days_antibiotics_therapy", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

rabu.bac.iap <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Total_days_antibiotics_therapy", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

ggarrange(rabu.bac.iap,rabu.vir.iap, common.legend = TRUE, legend = "bottom")

11.5 Antibiotics at birth

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.tdab <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Antibiotics_at_birth", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

rabu.bac.tdab <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Antibiotics_at_birth", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:3]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

ggarrange(rabu.bac.tdab,rabu.vir.tdab, common.legend = TRUE, legend = "bottom")

11.6 Small for gestational age

#remotes::install_github("jstokholm/rabuplot")
library(rabuplot)

rabu.vir.sga <- rabuplot(
  PSV.host,
  violin = TRUE, predictor = "Small_for_gestational_age", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:2]) +
  ggtitle("Viral hosts") +
  scale_x_discrete(limits = rev)

rabu.bac.sga <- rabuplot(
  PSB,
  violin = TRUE, predictor = "Small_for_gestational_age", p_adjust = TRUE, N_taxa = 10, colors = col_fil[1:2]) +
  #theme(legend.position = "none") +
  ggtitle("Bacteria") +
  scale_x_discrete(limits = rev)

ggarrange(rabu.bac.sga,rabu.vir.sga, common.legend = TRUE, legend = "bottom")

12 Differential abundance - Deseq2

12.1 Days

12.1.1 Bacteria

12.1.2 Virome

12.1.3 Virome hosts

12.2 BPD

12.2.1 Bacteria

12.2.2 Virome

12.2.3 Viral hosts

12.3 Small for gestational age

12.3.1 Bacteria

12.3.2 Virome

12.3.3 Viral hosts

13 Cross correlation

13.1 rCCA

13.1.1 Conversion and screeplot

#Subset samples present in both
vir_ids <- PSV@sam_data$description %>% str_split_i("_",-1)
PSV@sam_data$id <- vir_ids

PSB.CSS.both <- subset_samples(PSB,description %in% vir_ids)
#PSV.CSS.both <- subset_samples(PSV.CSS,vir_ids %in% PSB.CSS.both@sam_data$description)
PSV.CSS.host.both <- subset_samples(PSV.CSS.host,vir_ids %in% PSB.CSS.both@sam_data$description)

##Transpose and convert to dataframes
botu <- data.frame(PSB.CSS.both@otu_table)
#Select last non-empty taxonomic rank
tax.bac <- as.data.frame(PSB.CSS.both@tax_table) %>% dplyr::mutate_all(~na_if(., '')) %>% apply(.,1, function(x) tail(stats::na.omit(x), 1))
rownames(botu) <- paste(rownames(botu),tax.bac, sep = "-")

votu <- data.frame(PSV.CSS.host.both@otu_table)
#Select last non-empty taxonomic rank
tax.vir <- as.data.frame(PSV.CSS.host.both@tax_table) %>% dplyr::mutate_all(~na_if(., '')) %>% apply(.,1, function(x) tail(stats::na.omit(x), 1))
rownames(votu) <- paste(rownames(votu),tax.vir, sep = "-")


  ###########Convert to matrix and transpose - Maybe skip?
t(as.matrix.data.frame(botu))->tmbacteria
  #Save OTU Taxonomy from header
tmcolnames <- colnames(tmbacteria)
  #Remove first line
rownames(tmbacteria) <- NULL
  #Apply as.numeric to matrix(will remove header as it is in characters)
nbacteria <- mapply(tmbacteria, FUN=as.numeric)
  #Convert back into matrix and add colnames
mbacteria <- matrix(nbacteria, ncol=NCOL(tmbacteria), nrow=NROW(tmbacteria))
colnames(mbacteria) <- tmcolnames

  
    ##Repeat for viral OTU table
t(as.matrix.data.frame(votu))->tmvirome
tmcolnames <- colnames(tmvirome)
rownames(tmvirome) <- NULL
nvirome <- mapply(tmvirome, FUN=as.numeric)
mvirome <- matrix(nvirome, ncol=NCOL(tmvirome), nrow=NROW(tmvirome))
colnames(mvirome) <- tmcolnames
    ##Alternatively use: data.matrix(frame, rownames.force = NA)

    ##Converting and removing near-zero variance values
bac <- mbacteria[,-c(1)]
nzv <- nearZeroVar(bac)
bac <- bac[, -nzv$Position]
           
           
vir <- mvirome[,-c(1)]
#nzv <- nearZeroVar(vir) # If any
#vir <- vir[, -nzv$Position]


  #PCA
library(vegan)
mod_bac <- rda(bac, scale = TRUE)
plot(mod_bac, scaling = 3)

screeplot(mod_bac)

mod_vir <- rda(vir, scale = TRUE)
plot(mod_vir, scaling = 3)

screeplot(mod_vir)

scaling = "shrinkage"
  #Apply scaling 
result_1comp <- rcc(bac, vir, method = scaling, ncomp = 1)
result_2comp <- rcc(bac, vir, method = scaling, ncomp = 2)
result_3comp <- rcc(bac, vir, method = scaling, ncomp = 3)
result_4comp <- rcc(bac, vir, method = scaling, ncomp = 4)

ggarrange(plot(result_1comp, type = "barplot"),
          plot(result_2comp, type = "barplot"),
          plot(result_3comp, type = "barplot"),
          plot(result_4comp, type = "barplot")
          ,nrow = 2)

## $`1`

## 
## $`2`

## 
## attr(,"class")
## [1] "list"      "ggarrange"

# produce a heat map of the cross correlation matrix
imgCor(bac, vir, sideColors = c("purple", "green")) 

13.1.2 Model tuning

# set grid search values for each regularisation parameter
grid1 <- seq(0.001, 0.2, length = 20) 
grid2 <- seq(0.001, 0.2, length = 20)

# optimise the regularisation parameter values
sl(cv.tune.rcc, tune.rcc(bac, vir, grid1 = grid1, grid2 = grid2, 
                                   validation = "loo"),overwrite = F)

# examine the results of CV tuning
plot(cv.tune.rcc)

cv.tune.rcc

## 
## Call:
##  tune.rcc(X = bac, Y = vir, grid1 = grid1, grid2 = grid2, validation = "loo") 
## 
##   lambda1 =  0.1057368 , see object$opt.lambda1
##   lambda2 =  0.01147368 ,  see object$opt.lambda2
##  CV-score =  0.2003135 , see object$opt.score

#manually set lambda scores
#opt.l1 <- 1.7
#opt.l2 <- 0.003
opt.l1 <- cv.tune.rcc$opt.lambda1
opt.l2 <- cv.tune.rcc$opt.lambda2

13.1.3 Model and component selection

shrink.rcc <- rcc(bac, vir, method = "shrinkage")

CV.rcc <- rcc(bac, vir, method = "ridge", 
                         lambda1 = opt.l1,
                         lambda2 = opt.l2) 

# barplot of cross validation method rCCA canonical correlations
plot(CV.rcc, type = "barplot", main = "Cross Validation") 

# barplot of shrinkage method rCCA canonical correlations
plot(shrink.rcc, type = "barplot", main = "Shrinkage")

13.1.4 Run model

#Perfrom optimised CV rCCA

result.CV <- rcc(bac, vir, method = "ridge",
                         lambda1 = opt.l1,
                         lambda2 = opt.l2,
              ncomp=3)

result.raw <- rcc(bac, vir, method = "shrinkage",ncomp=3)

#Remove OTU IDs to only show taxonomy
result <- result.raw

result$names$colnames$X <- gsub(".*\\-", "", result$names$colnames$X) %>% make.unique(sep="_")
result$names$colnames$Y <- gsub(".*\\-", "", result$names$colnames$Y) %>% make.unique(sep="_")


# colnames(result$X) <- gsub(".*\\-", "", colnames(result$X))
# colnames(result$Y) <- gsub(".*\\-", "", colnames(result$Y))

# rownames(result$loadings$X) <- gsub(".*\\-", "", rownames(result$loadings$X))
# rownames(result$loadings$Y) <- gsub(".*\\-", "", rownames(result$loadings$Y))

13.1.5 Visualization

13.1.5.1 Heatmaps

#Defne pallette to only add colour to only high
##Generate new pallette
palette <- hcl.colors(200, palette = "Zissou 1", alpha = NULL, rev = FALSE, fixup = TRUE) #Seclect palette https://developer.r-project.org/Blog/public/2019/04/01/hcl-based-color-palettes-in-grdevices/
palette[61:140] <- "white" #Select which part of pallette to replace with white

# high4_2c <- cim(result, comp = 1:2, ylab = "Bacteria", xlab = "Virome", margins = c(12,18),cutoff = 0.4, color = palette)
# high45_2c <- cim(result, comp = 1:2, ylab = "Bacteria", xlab = "Virome", margins = c(12,18),cutoff = 0.45, color = palette)
#high5_3c <- cim(result, comp = 1:2, ylab = "Bacteria", xlab = "Virome", margins = c(12,18),cutoff = 0.5, color = palette)

#pdf(filename="high5_3c.pdf", res = 600, width = 4000, height = 4000)
#high5_3c
#dev.off()

######

high6_3c <- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(7,9),cutoff = 0.6, color = palette)

high5_3c <- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(7,9),cutoff = 0.5, color = palette)

high45_3c<- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(15,20),cutoff = 0.45, color = palette)

high4_3c <- cim(result, comp = 1:3, ylab = "Bacteria", xlab = "Virome host prediction", margins = c(15,20),cutoff = 0.4, color = palette)

#mat5 <- high5_2c$mat
mat6_3c <- high6_3c$mat
mat45_3c <- high45_3c$mat
mat4_3c <- high4_3c$mat
  
 #Output Matrix
write.table(mat45_3c, "Matrix_high_3-Component_0.45.txt", sep="\t")

#Make linear plots of all bac and vir 
bacmat6_3c <- bac[,rownames(mat6_3c)]
virmat6_3c <- vir[,colnames(mat6_3c)]
bacvirmat6_3c <- cbind(bacmat6_3c,virmat6_3c)
bac45_3c <- bac[,rownames(mat45_3c)]
vir45_3c <- vir[,colnames(mat45_3c)]
bacvirmat45_3c <- cbind(bac45_3c,vir45_3c)

#ggpairs
library(GGally)
#ggpairs(bacvir)
#ggpairs(bacvir45,columns=3:13)
#ggpairs(as.data.frame(bacvir45))
sl(pairs45,ggpairs(as.data.frame(bacvirmat45_3c)))
sl(pairs6,ggpairs(as.data.frame(bacvirmat6_3c)))

# pairs45 +
#   theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) 

#pairs6 +
#   theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5)) 

13.1.5.2 Network plots

13.1.5.2.1 Shrinkage

# Reformat for gephi

## Save blocks as polygon to determine shape in network

result$names$Polygon <- result$names$data
result$names$Polygon <- replace(result$names$data, c("bac","vir"),c(3,4)) %>% as.integer()

## Remove ubique markers to allow shared coloring in Gephi
result$names_simple$colnames$X <- gsub("*.\\_", "", result$names$colnames$X)
result$names_simple$colnames$Y <- gsub("*.\\_", "", result$names$colnames$Y)

#Show network

network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cex.node.name = 14,
        cutoff = 0.3,
        color.node = col_fil[1:2]) 

network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.4,
        color.node = col_fil[1:2]) 

network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.5,
        color.node = col_fil[1:2]) 

#Export network result
network.res <- network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.30,
        color.node = col_fil[1:2])

library(igraph)

#network.res 

write.graph(network.res$gR, file = "network.shrinkage0.3.gml", format = "gml")

13.1.5.2.2 CV

result.CV$names_simple$colnames$X <- gsub("*.\\_", "", result.CV$names$colnames$X)
result.CV$names_simple$colnames$Y <- gsub("*.\\_", "", result.CV$names$colnames$Y)

#Show network
network(result.CV, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.3,
        color.node = col_fil[1:2]) 

network(result.CV, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.35,
        color.node = col_fil[1:2]) 

# network(result.CV, comp = 1:3, interactive = FALSE,
#         lwd.edge = 2,
#         cutoff = 0.5,
#         color.node = col_fil[1:2]) 


library(igraph)
network.res <- network(result, comp = 1:3, interactive = FALSE,
        lwd.edge = 2,
        cutoff = 0.30,
        color.node = col_fil[1:2])

write.graph(network.res$gR, file = "network.CV0.3.gml", format = "gml")

13.2 PROCrustes comparison matrix correlation quantifications

13.2.1 Subset matrices

bac0 <- bac[which(PSB.CSS.both@sam_data$Time == "T0"),] %>% t()
vir0 <- vir[which(PSV.CSS.host.both@sam_data$Time == "T0"),] %>% t()

bac15 <- bac[which(PSB.CSS.both@sam_data$Time == "T1"),] %>% t()
vir15 <- vir[which(PSV.CSS.host.both@sam_data$Time == "T1"),] %>% t()

bac30 <- bac[which(PSB.CSS.both@sam_data$Time == "T2"),] %>% t()
vir30 <- vir[which(PSV.CSS.host.both@sam_data$Time == "T2"),] %>% t()

bac90 <- bac[which(PSB.CSS.both@sam_data$Time == "T3"),] %>% t()
vir90 <- vir[which(PSV.CSS.host.both@sam_data$Time == "T3"),] %>% t()

data <- list(bac0,bac15,bac30,bac90,vir0,vir15,vir30,vir90)

13.2.2 PCoAs

get_pcoa <- function(ps="PSB.CSS"){

otu <- data.frame(ps@otu_table) %>% t()
bray <- vegan::vegdist(otu, method = "bray")
add <-  !(ade4::is.euclid(bray))
pcoa <- cmdscale(bray, k = nrow(otu)-1, eig = TRUE, add = add)
#
return(pcoa)
}

pcoa.bac <- get_pcoa(PSB.CSS.both)
vegan::ordiplot(pcoa.bac, type = "text", main = "PCoA for Bacteria Data, Bray Distances")

## species scores not available

pcoa.vir <- get_pcoa(PSV.CSS.host.both)
vegan::ordiplot(pcoa.vir, type = "text", main = "PCoA for Viral Data, Bray Distances")

## species scores not available

13.2.3 Procrustes

13.2.3.1 All samples

#All time points
pcoa.bac <- PSB.CSS.both %>% get_pcoa()
pcoa.vir <- PSV.CSS.host.both %>% get_pcoa()

sl(procrustes.test.all,vegan::protest(X = pcoa.bac, Y = pcoa.vir, scores = "sites", symmetric = T, permutations = 999999), overwrite = F)

procrustes.test.all

## 
## Call:
## vegan::protest(X = pcoa.bac, Y = pcoa.vir, scores = "sites",      permutations = 999999, symmetric = T) 
## 
## Procrustes Sum of Squares (m12 squared):        0.4219 
## Correlation in a symmetric Procrustes rotation: 0.7603 
## Significance:  0.029147 
## 
## Permutation: free
## Number of permutations: 999999

13.2.3.2 Per time point

#Individual time points
pcoa.bac0 <- PSB.CSS.both %>% subset_samples(Days == 1) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)
pcoa.bac15 <- PSB.CSS.both %>% subset_samples(Days == 15) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)
pcoa.bac30 <- PSB.CSS.both %>% subset_samples(Days == 30) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)
pcoa.bac90 <- PSB.CSS.both %>% subset_samples(Days == 90) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)

pcoa.vir0 <- PSV.CSS.host.both %>% subset_samples(Days == 1) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)
pcoa.vir15 <- PSV.CSS.host.both %>% subset_samples(Days == 15) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)
pcoa.vir30 <- PSV.CSS.host.both %>% subset_samples(Days == 30) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)
pcoa.vir90 <- PSV.CSS.host.both %>% subset_samples(Days == 90) %>% otu_table() %>% as.matrix() %>% as.data.frame() %>% filter(rowSums(.) != 0) %>% t() %>% vegan::vegdist(method="bray") %>% vegan::wcmdscale(eig = TRUE)

sl(procrustes.test0,vegan::protest(X = pcoa.bac0, Y = pcoa.vir0, scores = "sites", symmetric = T, permutations = 999999), overwrite = F)
sl(procrustes.test15,vegan::protest(X = pcoa.bac15, Y = pcoa.vir15, scores = "sites", symmetric = T, permutations = 999999), overwrite = F)
sl(procrustes.test30,vegan::protest(X = pcoa.bac30, Y = pcoa.vir30, scores = "sites", symmetric = T, permutations = 999999), overwrite = F)
sl(procrustes.test90,vegan::protest(X = pcoa.bac90, Y = pcoa.vir90, scores = "sites", symmetric = T, permutations = 999999), overwrite = F)


procrustes.test0

## 
## Call:
## vegan::protest(X = pcoa.bac0, Y = pcoa.vir0, scores = "sites",      permutations = 999999, symmetric = T) 
## 
## Procrustes Sum of Squares (m12 squared):        0.006186 
## Correlation in a symmetric Procrustes rotation: 0.9969 
## Significance:  0.33333 
## 
## Permutation: free
## Number of permutations: 5

procrustes.test15

## 
## Call:
## vegan::protest(X = pcoa.bac15, Y = pcoa.vir15, scores = "sites",      permutations = 999999, symmetric = T) 
## 
## Procrustes Sum of Squares (m12 squared):        0.3786 
## Correlation in a symmetric Procrustes rotation: 0.7883 
## Significance:  0.34134 
## 
## Permutation: free
## Number of permutations: 999999

procrustes.test30

## 
## Call:
## vegan::protest(X = pcoa.bac30, Y = pcoa.vir30, scores = "sites",      permutations = 999999, symmetric = T) 
## 
## Procrustes Sum of Squares (m12 squared):        0.4974 
## Correlation in a symmetric Procrustes rotation: 0.7089 
## Significance:  0.79621 
## 
## Permutation: free
## Number of permutations: 999999

procrustes.test90

## 
## Call:
## vegan::protest(X = pcoa.bac90, Y = pcoa.vir90, scores = "sites",      permutations = 999999, symmetric = T) 
## 
## Procrustes Sum of Squares (m12 squared):        0.3262 
## Correlation in a symmetric Procrustes rotation: 0.8209 
## Significance:  0.097494 
## 
## Permutation: free
## Number of permutations: 999999

#Combine  procrustes ressults in dataframe
test0 <- c(procrustes.test0$ss,procrustes.test0$scale,procrustes.test0$signif)
test15 <- c(procrustes.test15$ss,procrustes.test15$scale,procrustes.test15$signif)
test30 <- c(procrustes.test30$ss,procrustes.test30$scale,procrustes.test30$signif)
test90 <- c(procrustes.test90$ss,procrustes.test90$scale,procrustes.test90$signif)

m <- matrix(nrow = 0, ncol =3)
colnames(m) <- c("Sum of Squares", "Procrustes Correlation", "Pr.F")

m <- rbind(m,test0,test15,test30,test90)
m

##        Sum of Squares Procrustes Correlation      Pr.F
## test0     0.006185679              0.9969024 0.3333333
## test15    0.378600098              0.7882892 0.3413410
## test30    0.497403537              0.7089404 0.7962050
## test90    0.326183086              0.8208635 0.0974940

rownames(m) <- c("Day 0", "Day 15", "Day 30", "Day 90")

#Prepare for plotting

dat <- m %>%
  as.data.frame() %>%
  tibble::add_column(Variable = rownames(m),.name_repair = "minimal") 

#Add adjusted p-values
dat$Pr.F.adj <- p.adjust(dat$Pr.F, method = 'holm')

#Show table
#knitr::kable(anova.dat)

##Fancy plot

melt <- tidyr::pivot_longer(dat,cols=c("Procrustes Correlation"), names_to='measure', values_to="value") %>%
  add_significance(p.col = "Pr.F", output.col = "PrF.sym",
                   cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05, 0.1,1),
                   symbols = c("****", "***", "**", "*", "\u02D9", "")
                   ) #%>%
  #arrange(value)

procrustes.bar <- ggbarplot(melt,
                                         x="Variable",
                                         y="value",
                                         #color = "Variable",
                                         orientation = c("horizontal"),
                                         color = NA,
                                         fill = "measure",
                                         position=position_dodge(0.7),
                                         label = melt$PrF.sym,
                                         lab.size = 9,
                                         lab.vjust = 0.7,
                                         lab.hjust = -0.3
                                         
) +
  #stat_pvalue_manual(
  #  adonis.melt,  label = "PrF.sym", tip.length = 0.01
  #) +
  #geom_hline(yintercept = 0.05) +
  ggtitle("") +
  xlab("") +
  ylab("Procrustes test significance") +
  scale_fill_jco() +
  scale_x_discrete(limits = rev)

procrustes.bar

14 Figures

14.1 Figure 1

14.2 Figure S1

14.3 Figure 2

14.4 Figure S2

14.5 Figure S3

14.6 Figure 3

14.7 Figure S4

14.8 Figure S4 - CD?

15 Session info

sessionInfo()

## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Danish_Denmark.utf8  LC_CTYPE=Danish_Denmark.utf8   
## [3] LC_MONETARY=Danish_Denmark.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=Danish_Denmark.utf8    
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] igraph_1.5.1                GGally_2.2.0               
##  [3] rabuplot_0.0.1.07           UpSetR_1.4.0               
##  [5] VennDiagram_1.7.3           futile.logger_1.4.3        
##  [7] reshape2_1.4.4              metagenomeSeq_1.38.0       
##  [9] RColorBrewer_1.1-3          glmnet_4.1-8               
## [11] Matrix_1.6-1.1              limma_3.52.4               
## [13] readxl_1.4.3                phyloseq_1.40.0            
## [15] mixOmics_6.20.0             MASS_7.3-60                
## [17] usdm_2.1-7                  terra_1.7-65               
## [19] directlabels_2023.8.25      vegan_2.6-4                
## [21] lattice_0.21-8              permute_0.9-7              
## [23] viridis_0.6.4               viridisLite_0.4.2          
## [25] DESeq2_1.36.0               SummarizedExperiment_1.26.1
## [27] Biobase_2.56.0              MatrixGenerics_1.8.1       
## [29] matrixStats_1.0.0           GenomicRanges_1.48.0       
## [31] GenomeInfoDb_1.32.4         IRanges_2.30.1             
## [33] S4Vectors_0.34.0            BiocGenerics_0.42.0        
## [35] data.table_1.14.8           RVAideMemoire_0.9-83-7     
## [37] cowplot_1.1.2               plyr_1.8.8                 
## [39] ampvis2_2.8.6               rstatix_0.7.2              
## [41] ggsci_3.0.0                 ggpubr_0.6.0               
## [43] lubridate_1.9.2             forcats_1.0.0              
## [45] stringr_1.5.1               dplyr_1.1.3                
## [47] purrr_1.0.2                 readr_2.1.4                
## [49] tidyr_1.3.0                 tibble_3.2.1               
## [51] ggplot2_3.4.4               tidyverse_2.0.0            
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3          ragg_1.2.5              bit64_4.0.5            
##   [4] knitr_1.45              multcomp_1.4-25         DelayedArray_0.22.0    
##   [7] KEGGREST_1.36.3         MicrobiotaProcess_1.8.2 RCurl_1.98-1.12        
##  [10] generics_0.1.3          microbiome_1.18.0       lambda.r_1.2.4         
##  [13] TH.data_1.1-2           usethis_2.2.2           RSQLite_2.3.1          
##  [16] bit_4.0.5               tzdb_0.4.0              ggtreeExtra_1.6.1      
##  [19] httpuv_1.6.11           xfun_0.40               hms_1.1.3              
##  [22] jquerylib_0.1.4         evaluate_0.23           promises_1.2.1         
##  [25] fansi_1.0.4             caTools_1.18.2          DBI_1.2.0              
##  [28] geneplotter_1.74.0      htmlwidgets_1.6.4       rARPACK_0.11-0         
##  [31] ellipsis_0.3.2          RSpectra_0.16-1         ggnewscale_0.4.9       
##  [34] backports_1.4.1         annotate_1.74.0         libcoin_1.0-10         
##  [37] vctrs_0.6.3             remotes_2.4.2.1         here_1.0.1             
##  [40] abind_1.4-5             httr2_1.0.0             cachem_1.0.8           
##  [43] withr_2.5.2             treeio_1.20.2           cluster_2.1.4          
##  [46] ape_5.7-1               lazyeval_0.2.2          crayon_1.5.2           
##  [49] ellipse_0.5.0           genefilter_1.78.0       pkgconfig_2.0.3        
##  [52] labeling_0.4.3          nlme_3.1-163            pkgload_1.3.3          
##  [55] gitcreds_0.1.2          devtools_2.4.5          rlang_1.1.1            
##  [58] lifecycle_1.0.4         miniUI_0.1.1.1          sandwich_3.1-0         
##  [61] cellranger_1.1.0        rprojroot_2.0.4         aplot_0.2.2            
##  [64] raster_3.6-26           carData_3.0-5           boot_1.3-28.1          
##  [67] Rhdf5lib_1.18.2         zoo_1.8-12              png_0.1-8              
##  [70] bitops_1.0-7            KernSmooth_2.23-22      rhdf5filters_1.8.0     
##  [73] Biostrings_2.64.1       blob_1.2.4              shape_1.4.6            
##  [76] coin_1.4-3              gridGraphics_0.5-1      ggsignif_0.6.4         
##  [79] scales_1.3.0            memoise_2.0.1           magrittr_2.0.3         
##  [82] gplots_3.1.3            zlibbioc_1.42.0         compiler_4.2.1         
##  [85] lme4_1.1-34             ggstar_1.0.4            cli_3.6.1              
##  [88] ade4_1.7-22             XVector_0.36.0          urlchecker_1.0.1       
##  [91] lmerTest_3.1-3          patchwork_1.1.3         formatR_1.14           
##  [94] mgcv_1.9-0              tidyselect_1.2.0        stringi_1.7.12         
##  [97] textshaping_0.3.6       highr_0.10              yaml_2.3.7             
## [100] locfit_1.5-9.8          ggrepel_0.9.3           sass_0.4.8             
## [103] ggstats_0.5.1           tools_4.2.1             timechange_0.2.0       
## [106] parallel_4.2.1          rstudioapi_0.15.0       foreach_1.5.2          
## [109] gridExtra_2.3           Rtsne_0.16              farver_2.1.1           
## [112] digest_0.6.33           shiny_1.8.0             quadprog_1.5-8         
## [115] Rcpp_1.0.11             car_3.1-2               broom_1.0.5            
## [118] later_1.3.1             httr_1.4.7              AnnotationDbi_1.58.0   
## [121] Wrench_1.14.0           colorspace_2.1-0        XML_3.99-0.16          
## [124] fs_1.6.3                splines_4.2.1           yulab.utils_0.1.2      
## [127] tidytree_0.4.6          gh_1.4.0                sp_2.1-2               
## [130] multtest_2.52.0         ggplotify_0.1.2         systemfonts_1.0.4      
## [133] plotly_4.10.3           sessioninfo_1.2.2       xtable_1.8-4           
## [136] nloptr_2.0.3            jsonlite_1.8.7          ggtree_3.4.4           
## [139] futile.options_1.0.1    corpcor_1.6.10          modeltools_0.2-23      
## [142] ggfun_0.1.3             R6_2.5.1                profvis_0.3.8          
## [145] metagMisc_0.0.4         pillar_1.9.0            htmltools_0.5.7        
## [148] mime_0.12               minqa_1.2.6             glue_1.6.2             
## [151] fastmap_1.1.1           BiocParallel_1.30.4     codetools_0.2-19       
## [154] pkgbuild_1.4.3          mvtnorm_1.2-3           utf8_1.2.3             
## [157] bslib_0.6.1             numDeriv_2016.8-1.1     curl_5.2.0             
## [160] gtools_3.9.4            survival_3.5-7          rmarkdown_2.25         
## [163] biomformat_1.24.0       munsell_0.5.0           rhdf5_2.40.0           
## [166] GenomeInfoDbData_1.2.8  iterators_1.0.14        gtable_0.3.4